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ABSTRACT 



The objective is to obtain and compare interpolation 
methods and contour plot techniques , based on data from 
sites whose locations are irregular or scattered. Data 
from the 1980 Fugro report ^Ref. lj on MX valley soil 
samples serves as a test bed. 

Two groups of data, seismic p-wave velocities and sur- 
face soil depth were studied using six different interpola- 
tion methods and three different contour plot techniques. 
Comparison of the interpolation methods and of the contour 
plot techniques are made. 
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INTRODUCTION 



I . 

Detailed investigations of seismic p-wave velocities and 
depth, of the surface soil layer were performed in the Ralston 
Valley, Nevada, by Fugro , Inc.[_Ref. lj , for the U.S. Army 
Engineer Waterways Experiment Station (W.E.S.) as part of the 
preliminary work on the MX missile system deployment. Six- 
teen locations were chosen in the Ralston Valley by W.E.S. 
for their study organized into four sets of sites (coded RA, 
RB,RC,RD) located in each of the four predominant surficial 
soil types (coded 5I,5Y,U,4U), and an additional site marked 
RU2 identified in the Valley. 

The data source for identified sites [Ref. 2j gives the 
(X,Y) coordinates (one mile to the inch) of the boring 
locations relative to the map's lower left grid mark, the 
seismic p-wave velocities (MPS - meters per second) of the 
surface layer, and the depth of that layer (meters) for the 
2B line; see Appendix A. The site identifiers and locations 
for the Ralston Valley sampling plan can be found in Appendix 
B. The immediate concern is to develop contours describing 
the depth of the top layer, and of the compression wave speed 
in the top layer, in the entire Valley. 

Based on these limited data, the mathematical problem is 
that of constructing a smooth bivariate function F(X,Y) which 

) = F^ for K = 1,...,N. 



takes on specified values F(X^,Y^ 
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Franke |[Ref. 3,4,5, 6j has made a rather extensive general 
study of the modern interpolation methods. Most of the methods 
require the user to specify a parameter (or several) , but 
these parameters can be determined by computational experi- 
ment ^Ref . 4 J . The interpolation methods selected for 
specific inclusion in the current work are a subset of these 
methods, see Table 1. Furthermore, contour plotting, specific- 
ally non-IMSL subroutine CONTUR, subroutine CONISD and sub- 
routine PLT3D1 , were exercised for comparison. Basic FORTRAN 
programs for the first five interpolation methods were obtained 
from Professor R. Franke ([Ref. 3, 4,5,6 J. 

The various interpolation methods were used to supply 
function values at a gridwork of reference points. In order 
to compare the methods, the mean and standard deviation of 
the values were produced at each grid point. Also it is 
important to identify which methods generates sharply different 
interpolation values compared to the other methods. For this 
purpose, the same FORTRAN program generates the absolute dif- 
ferences from the mean values for each method. 

FORTRAN programs for all of the methods can be obtained 
by writing the advisor. 

In the next chapter we will describe a number of inter- 
polation methods. Chapter 3 describes employment of the 
Ralston Valley data as a way of comparing these methods. 

Chapter 4 discusses contour plot techniques, and the con- 
clusions from using these procedures are given in Chapter 5. 
Graphical results are presented in a series of appendices. 
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TABLE 1 



INTERPOLATION METHODS 

1) Inverse Distance Weighted Methods 

A) Shepard's Method 

B) Modified Shepard's Method 

C) Modified Linear Shepard's Method 

D) Modified Shepard's Method Boolean Sum Plane 

E) Modified McLain Method 

F) Quadratic Shepard's Method 

*G) Modified Quadratic Shepard's Method 
H) McLain Method 

2) Franke's Methods (Rectangle based blending methods) 

A) Franke's Method (Mode One) 

B) Franke's Method (Mode Three) 

*C) Franke's Method (Thin plate local functions) 

3) Triangle Based Blending Methods 

A) Nielson-Franke Linear Triangle Method 
*B ) Nielson-Franke Quadratic Triangle Method 

4) Finite Element Based Methods 

A) Akima's Method 

B) Akima's Method - Modification Two 

C) Akima's Method - Modification One 

D) Akima's Method - Modification Three 
*E) Nielson's Minimum Norm Network 

F) Lawson's Method 

5) Foley's Method 

A) Generalized Newton Interpolant 

B) TF Delta Sum Berstein Interpolant 

C) Iterated Delta Sums: TF Delta Sum Bicubic Spline 

*D) Iterated Delta Sums: A Shepard's Method Delta Sum 

Bic. Sp. 

6) Nodal Basis Function Type Methods 
A) Rotated Gaussian 

*B) Hardy's Multiquadric 

C) Hardy's Reciprocal Multiquadric 

D) Duchon's Radial Cubic Method 

E) Duchon's Thin Plate Function 

F) Rotated B-Spline 



* Methods chosen for comparison study. 
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II. INTERPOLATION METHODS 



The problem being addressed here is that of constructing 
a smooth bivariate function, F(X,Y) , which takes on specified 
values F(X^,Y k ) = F k , k=l,...,N based on point data from 
scattered, irregular locations. The bivariate function, 

F(X,Y) is smooth, (i.e., it has at least continuous first 
partial derivatives) , and the points (X^/Y^) are 'scattered', 
that is, not necessarily conforming to some regular pattern. 

Some interpolation methods, called 'global', are sensitive 
to all the data points : the addition or deletion of a data 

point, or a change of one of the coordinates of a data point, 
will propogate throughout the entire domain of definition. 

Other interpolation methods are called 'local', that is, insen- 
sitive to changes in the data points which are sufficiently 
distant from the location at hand. 

Franke, R. £Ref. 3, 4, 5, 6 J has made an extensive general 
study of new interpolation methods. He produced a classifi- 
cation of the interpolation methods ^Ref. 4j, as shown in 
Table 1. Under his six main types, there are a total of 29 
different methods. All the methods appearing in the same 
group are either local or global and use the same basic idea; 
several are modifications of the previously described ones. 

Six methods, one from each of the six main types, are 
described below. 
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A. MODIFIED QUADRATIC SHEPARD'S METHOD 



The modified quadratic Shepard's method was derived 
from the inverse distance weighted method. All methods of 
this type which are considered may be viewed as generaliza- 
tions of Shepard's method. The basic Shepard's method is 

N N 

F (X , Y) = Z W. (X ,Y) . Q, (X ,Y) / Z W. (X,Y) 
k=l * k=l * 



-y 

where W^(X,Y) = d^ . Typically y = 2, although other values 

may be used Here d 2 = ( (X-K^) 2 + (Y-Y^) 2 ). 

In the modified quadratic Shepard's method the weights 
for obtaining the nodal functions (quadratics) are taken as 

2 

and R = 1/2 I/n/N.'D 
q * 4 

where D is the diameter of the point set, and a nominal value 

for N , determined by computational experiment £Ref . 4 J[, is 
4 

N =18. The diameter of the point set is defined 
4 



W (X,Y) = 
x 



(R q - V 



(R„ 



V 



2 

D = MAX 
I,J 



2 2 
(X. - X.) + (Y. - Y.) 

ID 1 D 



The nodal functions Q, (X,Y) are defined as 

k 
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Q (X,Y) = F + a (X-X ) + a (Y-Y, ) + a (X-X ) 



k k2 k 



k3 k 



k 4 k 



+ a (X-X ) • (Y-Y ) + a, (Y-Y ) . 

k5 k k kg k 



and for the coefficients one must solve the least squares 
problem. 



MIN 



W * [ F k + **2 ( VV + WW 



a^ , j = 2 , . . . , 6 i£k 



+ a k4 (X i* X k ) + a k5 ( VV- (Y i-V 



+ a k6 (Y i- Y k> 2 ' F i] 



Complete details are given in[]Ref. 4 , 5 J . 

B. LOCAL THIN PLATE SPLINES METHOD 

This is known as Franke's method ^Ref. 4 J. The idea is 
to represent the interpolation function as 



F (X,Y) = Z W, . (X,Y) -Q. . (X ,Y) 



1/1 



1>1 






where the (X,Y) are local weight functions j^Ref. 6j. Here 
we choose the W so that W. . = 1, and the local approximations 






i,l 



Q^j (X,Y) interpolate at points where the corresponding weight 
function is non-zero. The weight functions are products of 
piecewise hermite cubics^Ref. sj. 
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The local approximations for this method are taken to be 



the thin plate splines and have the form 



Qi , (X,Y) 
“'vJ 



2 

Z A v *d, -LOG d.+a + b'X + C-Y 
k£l k k 



where dj“ = 



((X-X k ) + (Y-Y k ) ) 



And I is the set of indices k for which (X^,Y^,F k ) 
to be interpolated by Q i j(X,Y) , and is defined by 



is a point 



I 



K : Q is to take the value F^ at 




The coefficients A k , and a, b and c are determined by a certain 
linear system of equations, (see [^Ref. 6]) . 



C. N I ELS ON- F RANKE QUADRATIC TRIANGLE METHOD 

All methods listed under triangle based blending methods 
[ Ref. 4 J are conceptually similar to the inverse distance 
weighted methods. The interpolation function for the quad- 
ratic triangle method is 

F (X , Y) = W i (X / Y) ‘Q ± (X / Y) + W- (X,Y) 'Q- (X,Y) + W R (X, Y) -Q k (X,Y) 

The first step in the present method is to partition the plane 
into triangles by connecting neighboring data points based 
upon the min-max angle criterion as described by Franke in 
£*Ref. 5j. The quadratic triangle method uses the inverse 
distance weighted quadratic Q^(X,Y) which is also used in the 
modified quadratic Shepard's method, see section 2. A. 
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A significant difference from the inverse distance weighted 
method is that the weight functions for the triangle based 
blending method are based on a triangulation of the convex 



tive to the triangulation technique. If a triangulation is in 
existence for other purposes, it can be used. 

In short, the method is to define the nodal functions 
/k=l , . . . ,N , as in the modified quadratic Shepard's method, 
form a triangulation of the points , determine the weight 
functions ^Ref. sj and compute interpolant F(X,Y). 

D. NIELSON'S MINIMUM NORM NETWORK METHOD 

Nielson's minimum norm network method ^Ref. 4,5_], can 

be found under the finite element based methods in Table 1. 

Like the Nielson-Franke quadratic triangle method, this global 

method is based upon a tri angulation . 

The method consists of three separate steps; triangulation 

curve network, and blending J^Ref. 7_]. The points 

(X ,Y, ) ,k=l , . . . ,N are used as the vertices of a triangulation 
k K 

of the convex hull of these points. Alternatively, if a tri- 
aggulation already exists, it can be used. The curve network 
step involves the solution of a certain minimum norm problem 
and requires first partial derivatives in its discretized 
form. These are obtained by assuming a cubic variation along 
each edge in the triangulation and minimizing the integral of 
the second derivative squared ^Ref. 7_J. 



hull of the point set 




The method is not sensi- 
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This method does not provide extrapolation outside the con- 
vex hull, but the triangular blending method can be used to 
extend the curve network to the entire hull. In a triangular 
exterior region, the function is taken to be the linear 
function determined by the value and slopes at the vertex 
[Ref. ?]. 



E. GENERALIZED NEWTON DELTA SUM BICUBIC SPLINE METHOD 

Foley's methods [Ref. 4 J involve several ideas. The 
interpolant is taken as either the generalized Newton inter- 
polant, or a form of Shepard's method. The use of generalized 
Newton type interpolant is involved in them prominently, 
because the best performance for various test functions is 
provided by the interated delta sum methods using the generalized 
Newton interpolant with natural bicubic splines. 

The generalized Newton interpolant takes the form 

N 

T (X,Y) = I a, *W k (X,Y) r 
N k=l k * 



where 



a k 



f 

k 



T 

k-1 



(x ,y ) 
k k 



vw 



and W ( X , Y ) has the property W (X. ,Y.) = 0 i=l ,2 , . . . ,k-l . 
k k 1 1 

This function is dependent upon the order of the points, and 
so Foley's scheme involves an ordering process [Ref. 4_]. 
After constructing a bicubic spline interpolant for the grid 
points, one adds a generalized Newton interpolant for the 
difference between data and the bicubic splines, obtaining 
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an interpolant . This process is termed a 'delta sum' by 
Foley (^Ref. 4J. 



F. HARDY'S MULTIQUADRIC METHOD 

Hardy's multiquadric method [Ref. 8,9] can be thought 
of as being under the global basis function type. A multi- 
quadric surface can be represented by 



N 

Z 

j = l 



2 2 
(X.-X. ) + (Y -Y ) + C 

D 1 j i 



1/2 






i=l , . . . ,N 



The given problem data provide a set of cartesian coordinates 
on the surface ranging from X^,Y^,Z^ to X^/Y^/Z^; the quad- 
ratic term coefficients C-^ to are unknown. The coordinates 
of the N data points can be substituted into the definition 
of the surface and only the coefficients C. are left as 



D 

unknown. A problem here is that of solving a system of N 
linear equations with N unknowns . A nominal value for C was 
determined by computational experiment and the value C = 0 
was chosen. 



If a vector of unknowns is to be 
used, define 






then each known element becomes 
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1/2 



(X.-X. ) 
3 i 



(Y.-Y. ) 
3 i 



= a- 



ID 



and the matrix £ a j_j] = A, which is (NXN) 
Upon 

defining which t^il = 



coefficient matrix. 




i 



I 



A multiquadric surface reduces to AxC = Z and, as usual, the 
solution is C = A~^Z. 

When the known constraints C. are substituted into the 

D 

multiquadric surface formula, we have the required equation 
of a surface, which fits the data points exactly and which 
provides logical interpolation at intermediate points. 
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III. COMPARISON OF INTERPOLATION METHODS USING 
RALSTON VALLEY DATA 



Contour construction for seismic p-wave velocities and 
depth data can be seen in Appendix C and in Appendix F , 
respectively, for the various interpolation methods and with 
the three different contour plot techniques. Each method 
appears in the order given in Table 1. Each method is dis- 
played using first the subroutine CONTUR, then the subroutine 
CON ISP , and lastly the subroutine PLT3D1 . Subroutine CONTUR 
cannot outline the Ralston Valley, but gives contour con- 
structions in rectangular form with the x-axis ranging from 
one to seventeen miles and the y-axis from one to twenty five 
miles. The subroutine CONISD can outline the Ralston Valley 
with boundary lines. All three dimensional pictures were 
generated by subroutine PLT3D1. 

Appendix D shows the mean and standard deviation of the 
various contour constructions for the seismic p-wave velo- 
cities. Numerical values for the mean and standard deviation 
at each grid point come from six different interpolation 
methods. Appendix G is similar to Appendix D, but it is for 
the depth of the surface soil layer. Convex contour lines 
are formed for the standard deviations values as can be seen 
in Appendices D and G. Small standard deviation values 
appear in data point region, and become larger as one moves 
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toward the edges. In other words, more accurate prediction 
occurs in data point region and less accurate prediction out- 
side this region. 

Deviations from mean values for all interpolation methods 
can be seen in Appendix E for seismic p-wave velocities and 
in Appendix H for depth values. The 1, 10, 50 and 100 contours 
show the absolute differences between the given method and 
mean values for seismic p-wave velocities. The 0.1, 0.5, 1.0 
contours for depth data display quite convex regions for all 
methods. Clearly, it can be seen that*, there are no signifi- 
cant differences between methods in terms of deviations from 
mean values for grid points especially in data region. This 
result shows that in the data points region all methods work 
with almost same accuracy and outside the region there are 
too great differences. There is no way to compare them out- 
side the region. 
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IV. CONTOUR PLOT TECHNIQUES 



The interpolation values provided by any of the methods 
can be seen in the computer outputs , but their interpretation 
is highly difficult without contour construction. For inter- 
pretation purposes , computer programs for drawing plots on 
an off-line plotter are available at the Naval Postgraduate 
School computer center. 

The subprograms CONTUR, CONISD and PLT3D1 serve the com- 
mon purpose but have different capabilities. 

All the subprograms have the same limitation in that they 
use a linear interpolation process to find the values of the 
points along contour segments. Also, there is a limit in the 
resolution of the off-line plotter, namely, .005 inch in both 
the X and Y directions for the versatec plotter. 

A. SUBROUTINE CONTUR 

Given a matrix of numerical values of Z = (X,Y) , the 
program CONTUR generates a contour graph on which one or 
more contour levels are drawn. This non-IMSL subroutine 
uses a two-dimensional array and turns out a two-dimensional 
off-line plot. Interior and exterior contour segments and 
labels (if requested) for exterior contour segments can be 
seen. Labeling of the interior segments represent local 
maxima. 
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Advantages : 

1) The subroutine has a minimum number of arguments (10) 
compared with other plot techniques. 

2) A rather big two-dimensional array can be used, up to 
1797 data points. 

3) It takes approximately one fiftieth of the CPU time that 
would be required by subroutine CONISD. 

Disadvantages : 

1) Irregularly spaced data cannot be contoured. 

2) Irregular boundary capability is not available, and there 
is no option for a 'cut-out' area. 

3) The width and the height of the contour graph must be 
specified as integers, and are not adjustable scale 
parameters . 

4) The scale values (if requested) one inch apart on the 
exterior frame of the contour graph increase from north 
to south. This is not oriented properly for our use, 
which is a lower left grid oriented system. 

5) Limited labeling. The x-axis and y-axis cannot be 
labeled. 

6) There is no option for smoothing the contour lines. 

3. SUBROUTINE CONISD 

This non-IMSL subroutine produces a drawn contour map on 

an off-line plotter for given irregularly spaced data points. 

Each data point is a triad of X,Y and Z values where 

Z = (X , Y) . There may be one or more 'cut-out areas' which 



are not contour. 
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Advantages : 

1) Irregularly spaced data can be drawn. 

2) Irregular boundary capability is available for one or 
more 'cut-out areas'. 

3) The scale for the x-axis and y-axis are defined as real 
numbers, i.e. , adjustable scale. 

4) X-axis and y-axis can be labeled. 

5) There is an option for smoothing the contour lines. 
Disadvantages : 

1) It takes much more CPU time than the other plot techniques. 
C. SUBROUTINE PLT3D1 

A subroutine PLT3D1 produces three-dimensional perspec- 
tive or isometric plots on an off-line plotter for given a 
two-dimensional matrix of numerical values of Z = (X,Y) . 
Labeling and scaling are not available for this reason, it 
can be used only for general purposes. 

A memorandum about this subprogram QRef . 12 J can be 
requested from W.R. Church computer center. This memorandum 
includes a listing of the computer program, the algorithm 
that is used by the program, and the axes orientation, rota- 
tion and projection are clarified by illustrations. 

Advantages : 

1) This routine gives a nice overview to given area by 
Z - (X ,Y) . 

2) The vantage point can be changed. 

3) The scale is adjustable. 



27 



Disadvantages : 

1) There is no option for 'cut-out areas'. 

2) There is no labeling. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



Seismic p-wave velocities or depth values for any grid 
point in Ralston Valley can be found from computer outputs or 
contour constructions for any particular interpolation method 
or for the mean of the six different methods. The standard 
deviations for each grid point are also available. The means 
and standard deviations may be useful for the comparison of 
methods . 

The computer outputs for standard deviations in Appendix 
E and H show that accuracy is high in the data point region 
and at some undecided distance from this region; but this 
distance cannot be measured. In the data point region all 
interpolation methods gives essentially the same result, 
regardless of whether they are global or local. 

Outside the data point region all methods generate rather 
different values and there is no way to decide which one is 
close to true values. If the point of interest is in the data 
point region any of the interpolation methods can be used. 

According to Appendix C and F; none of the methods can 
find a maximum or minimum point outside the region. Table 2 
gives approximate CPU times for all of the studied inter- 
polation methods. If minimal computer time is required for a 
given problem, the triangle based quadratic method is suitable. 
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CPU TIMES FOR STUDIED INTERPOLATION METHODS 
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APPENDIX A 



TABLE 3 

DATA SELECTED FOR TEST BED 



Side 


Designation 


Coordinates 
X Y 


P-Wave Vel. 
NPS 


Depth 

M. 


1 . 


RA5I 


9.25 


30.25 


561 


3.0 


2. 


RB5I 


5.69 


25.5 


366 


1.0 


3. 


RC5I 


4.75 


20.81 


279 


0.9 


4. 


RD5I 


5.75 


14.75 


347 


1.5 


5. 


RA5Y 


5.59 


26.25 


341 


1.5 


6. 


RB5Y 


11.75 


24.19 


366 


0. 8 


7. 


RC5Y 


7.06 


20.0 


323 


1.2 


8. 


RD5Y 


7.31 


10.31 


329 


0.8 


9. 


RAU 


7.31 


28. 19 


363 


1.2 


10. 


RBU 


9.25 


19.31 


335 


1.6 


11. 


RCU 


8.75 


11.25 


317 


1.4 


12. 


RDU 


17.44 


10.94 


378 


1.5 


13. 


RA4U 


10.12 


23.44 


341 


1.5 


14. 


RB4U 


12.25 


18.56 


372 


1.5 


15. 


RC4U 


11.44 


14.38 


372 


1.5 


16. 


RD4U 


15.44 


12.50 


381 


0.5 


17. 


RU2 


6.00 


27.19 


357 


1.2 
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APPENDIX B 



EXPLANATION 
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Figure 1 , 
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APPENDIX C 



P-WRVE VEL. (NPS) 
QUADRATIC SHEPARDS METHOD 
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(Using subroutine CONISD) 



Figure 3 
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Figure 4 
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Figure 5 
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LQCPL THIN PLRTE SPLINE M. 

(Using subroutine PLT3D1) 



Figure 7 
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Figure 9 
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(Using subroutine PLT3D1) 



Figure 10 
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Figure 11 
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P-WflVE VEL. (NPS) 
NIELSON’S MIN. NORM NETWORK M. 

(Using subroutine PLT3D1) 



Figure 13 
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Figure 14 
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Figure 15 
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Figure 19 
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Figure 23 . 
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